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We present a calculation of the fully differential cross section for Higgs boson production in the 
gluon fusion channel through next-to-next-to-leading order in perturbative QCD. We apply the 
method introduced in to compute double real emission corrections. Our calculation permits 
arbitrary cuts on the final state in the reaction hh —* H + X. ft can be easily extended to include 
decays of the Higgs boson into observable final states. In this Letter, we discuss the most important 
features of the calculation, and present some examples of physical applications that illustrate the 
range of observables that can be studied using our result. We compute the NNLO rapidity distri- 
bution of the Higgs boson, and also calculate the NNLO rapidity distribution with a veto on jet 
activity. 



PACS numbers: 

With Run II of the Tevatron producing data and the 
LHC set to begin operation in 2007, hadron colliders will 
soon play a major role in understanding the mechanism 
of electroweak symmetry breaking. Within the Stan- 
dard Model, this mechanism is linked to the Higgs bo- 
son, a scalar particle whose non-zero vacuum expectation 
value gives rise to the masses of all elementary particles. 
Finding the Higgs boson and analyzing its properties are 
therefore important parts of the high-energy physics pro- 
gram in the next decade. 

Hadron colliders, while offering substantial increases 
in energy over existing lepton machines, present several 
obstacles to performing precision physics studies. Non- 
perturbative QCD enters the calculation of cross sec- 
tions through both parton distribution functions (pdfs) 
of hadrons in the initial state and properties of jets in 
the final state. Perturbatively, the large value of the 
QCD coupling constant and the enhanced sensitivity to 
the factorization and renormalization scales force calcu- 
lations of many important processes to be extended to 
next-to-next-to-leading order (NNLO) in the QCD cou- 
pling constant. This reduces the unphysical sensitivity of 
the result to the renormalization and factorization scales; 
also, the additional partons model the structures of both 
colliding hadrons and final-state jets more accurately. 



There has been significant progress in the past several 
years in performing inclusive and semi-inclusive NNLO 
calculations y, |j, ||| . However, because of the cuts 
on the final states typical for the LHC and the Teva- 
tron, such results are of limited use. A fully differential, 
partonic Monte Carlo event generator is preferable. It 
is then guaranteed that, to a given order in the pertur- 
bative expansion, there is full control over hard emis- 
sions and the normalization of a given observable. How- 
ever, this approach can not be applied if the observable 
is dominated by regions of the phase space where soft 
and collinear effects are enhanced, and the resummation 
of certain higher order corrections may be necessary. In 
such situations, shower Monte Carlo event generators are 
used; they correctly describe the soft and collinear limits 
of the real emission. However, they can not reproduce 
the properties of hard radiation, and are unsuitable for 
precision studies. Ideally, we would combine perturbative 
results with shower Monte Carlos to gain the advantages 
of both approaches. This has been achieved at NLO 0], 
but has not been extended to higher orders. Construct- 
ing fully differential perturbative results is a first step 
toward this goal. 

Extending exclusive calculations to higher orders is not 
straightforward. While at NLO this problem has been 
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solved |8|, a similar solution at NNLO is not yet avail- 
able, despite significant effort 0- The primary obstacle 
is the double real radiation contribution to NNLO cross 
sections, which contains two additional partons in the fi- 
nal state. The singularity structure resulting from these 
emissions is significantly more complex than at NLO, and 
prevents the use of NLO techniques. 

We have recently suggested an alternative approach to 
the problem of real radiation at NNLO Q , which allows 
differential results to be obtained in an efficient, auto- 
mated fashion. We have tested our idea on the realistic 
example of e+e~ -> 2 jets at NNLO HE3 In this Let- 
ter we apply our method to calculate the fully differential 
cross section for Higgs boson production at hadron col- 
liders through NNLO in QCD. Since our result retains all 
kinematic information, it can be used to compute arbi- 
trary differential distributions, or to construct a partonic 
event generator accurate to NNLO. 

The dominant production mechanism for a light Higgs 
boson at hadron colliders is gluon fusion, gg — ► H, 
through a top quark loop. If the Higgs boson is suffi- 
ciently light, to/i < 2m t , its coupling to gluons can be 
described by a point-like vertex. Within this approx- 
imation, the NNLO corrections to the inclusive cross- 
section for Higgs hadroproduction have been considered 
recently in , where a detailed description of the setup 
can be found. The two-loop virtual corrections required 
for Higgs production are identical for differential distri- 
butions and the total cross section. The one-loop cor- 
rections to single real-radiation processes (gg — > H + g, 
qq — > H + g, qg — > H + g) can be computed easily with 
established methods The difficult contributions are 
the double real emission corrections that first appear at 
NNLO. 

We use the method first presented in to calculate 
these components. We combine an expansion in plus dis- 
tributions with sector decomposition to separate and 
extract their singularities. This requires a parameteriza- 
tion of the phase space in which the integration region 
is the unit hypercube. In principle, any mapping that 
accomplishes this is acceptable. In practice, finding a 
convenient parameterization that reduces the number of 
sector decompositions is important for the efficiency of 
the approach. 

For the double real emission corrections to Higgs pro- 
duction at NNLO, we must parameterize a 2 -> 3 par- 
ticle phase space, with one massive final-state particle. 
We consider here g(pi) + gip-z) — > H(h) + g(pz) + g{pi) 
as a prototypical partonic process. For a fixed energy of 
the partonic collision (jp\ + P2) 2 = s, the partonic phase 
space is described by four independent variables. We 
found it useful to employ several different parameteri- 
zations of the partonic phase space. We present here 
explicit formulae only for the parameterization which is 
the most convenient choice for the majority of diagrams. 
The scalar products s 34 = (p 3 + Pi) 2 , s ih = (pi ~Ph) 2 , 
and S2j = (p2 —Pj) 2 , where i = 1,2 and j = 3,4, take 
on simple forms in this parameterization, while the s\j 



contain overlapping singularities, and require sector de- 
composition when appearing in the denominator. In this 
parameterization, the 2^3 phase space becomes 
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where d = 4 — 2e is the space-time dimensionality, 
[dq] = dq d - 1 /{2q ), z = m 2 h /s, N = Q(d-2)Sl(d-3)(l- 
zf- 4e /2 4+2e , fl(d) = 27r d / 2 /T(d/2), and the independent 
scalar products are 

si h = -A 3 (l - z) [1 - Air(l - rt)/(r + «)] , 
S23 = -A 2 A 3 (l-z)[l + Ai(l-r<)/r/(r + i)], 
S34 = AiA 3 (l-A 3 )(l-z) 2 (l + u) 2 /^ P A, 
(l-A 3 )(l-z) 



Sl3 



K p r[l + Xi(l-rt)/r/(r + t)] 



A 1 +A 2 + 2(2X i -l)y/A 1 A 2 



(2) 



We have introduced the notations u = (z + A 3 (l — z))/(l — 
A 3 (l - z)), r = y/u, t = yfz, K P}Tn = (r ± i)(l ± 
rt), A 1 =X 1 {l-X 2 )(l+u) 2 , A 2 = A 2 (l - \i)r(K p - 
X\K m ). 

With the expressions given above, it is in principle 
straightforward to apply the method explained in to 
derive the finite, fully differential cross section for Higgs 
production at NNLO. However, some technical points 
warrant further discussion. A brute force application of 
the algorithm in 1] typically produces very large expres- 
sions. It is therefore important to organize the calcula- 
tion efficiently. 

• We found it important to identify all possible sym- 
metries of the process, and utilize the fact that 
many terms in the matrix elements are identical 
under a simple rotation of momentum labels, in or- 
der to reduce expression sizes. 

• We found it useful to introduce a specific phase- 
space parameterization only for the denominators 
of the matrix elements. We keep the numerators 
written in terms of invariant masses. We there- 
fore divide the expressions into universal denomi- 
nator structures (topologies), valid for any 2 — ► 3 
real emission correction with one massive final-state 
particle, and process-dependent numerators. This 
allows us to implement other processes in our code 
very simply. 

• A few topologies contain denominators that depend 
quadratically upon the A;. Since the full matrix el- 
ements do not contain quadratic singularities, there 
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must be numerator structures responsible for regu- 
lating this behavior. We therefore cannot trivially 
keep the numerator written in terms of invariant 
masses for these terms. It is usually simple to iden- 
tify the required numerator structure. An example 
found in Higgs production is the topology 

■A(s tj ) _ (S13S24 - S34) 2 , ,„x 

„2 „2 — „2 „2 >H s »iJ- V>) 

*23 6 14 *23*14 

The left-hand side of Eq. © naively has quadratic 
singularities, but, as shown, we can identify the 
required scalar product in the numerator that reg- 
ulates these. We can now write A(sij) using invari- 
ant masses, and introduce a specific parameteriza- 
tion for the remainder. 

• As in previously studied examples Q , the matrix el- 
ements contain "line" singularities that arise when 
a singularity is mapped to an edge of phase space. 
These are removed using an additional variable 
change A; — * A*, as discussed in 

• When we convolute the partonic cross sections with 
the pdfs to form the hadronic cross section, the 
variable z scales asz^ m\/ (xiXzShad), where Shad 
is the hadronic center-of-mass energy squared, and 
the Xi are the fractions of the hadronic momenta 
carried into the hard scattering process. The par- 
tonic cross sections are distributions in z, and are 
singular in the limit z — > 1. Is is clear from Eq.JQ) 
that the factor (1 — z)~ 4e which regulates this limit 
has been extracted, and therefore singularities in z 
and Xi can be treated identically. 

After the extraction of singularities, we can combine 
the double real emission corrections with the remaining 
contributions to the hard scattering cross section. This 
expression is then integrated numerically together with 
the pdfs to form the hadronic cross section. We perform 
the numerical integration with the version of VEGAS de- 
scribed in We use the MRST pdf sets, with mode=l. 
The numerical integration we perform is 6-dimensional; 
this includes the four independent partonic variables and 
the two Xi. The CPU times required for this calcula- 
tion depend strongly on both the kinematics and the re- 
quested precision. For Tevatron kinematics, reaching 1% 
precision on the inclusive cross section requires 0.5 hours, 
while obtaining 1% at the LHC requires about 1.5 hours. 
When constraints are imposed on the final state, the run 
times increase. 

We have performed several checks on our calculation. 
We cancel 1 /epoles numerically, as is typically done using 
this method 0,113. The singularities begin at 1/e , and 
involve all the different contributions to the final result: 
two-loop virtual corrections, ultraviolet renormalization, 
virtual corrections to NLO processes, double real emis- 
sion, and collinear subtractions for pdfs; the cancellation 
of the singularities is therefore a strong check on the cal- 
culation. We find that the singularities cancel to a high 
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FIG. 1: Bin-integrated rapidity distribution for LHC kine- 
matics. The scale [i is varied between m?,/2 < /i < 2mh- The 
LO, NLO, and NNLO distributions are shown. 



precision for both the inclusive cross section and all dis- 
tributions we have studied. For example, the e-expansion 
for the NNLO inclusive cross section at the LHC is 

a = (-0.05±2) x 10" 6 /e 4 + (-0-7±2) x 10" 3 /e 3 
+ (-0.3 ± 2) x 10" 3 /e 2 + (-3 ± 9) x 10~ 3 / e 
+ (44.9 ±0.4). (4) 

These results are in picobarns, and are obtained using 
rrih = 120 GeV and m t = 175 GeV; for this and other re- 
sults in this paper, higher order QCD corrections are ob- 
tained using the effective Lagrangian valid for nih < 2mt 
and then normalized to the exact tree-level cross sec- 
tion with its full mt-dependence. The renormalization 
and factorization scales have been set equal to the Higgs 
boson mass, /in — /if = fn>h- We have found a very 
good agreement between our result for the finite inclu- 
sive cross section for both Tevatron and LHC kinematics 
and several calculations of this quantity available in the 
literature Q. 

We now present several distributions that illustrate the 
range of observables that can be studied using our com- 
putation. The numerical precision of all results shown is 
1%. We first compute the bin-integrated rapidity distri- 
bution of the Higgs boson. We separate the entire rapid- 
ity range into twenty bins; because of the symmetry of 
the distribution under Y — > — Y, we only need to consider 
ten bins for Y > 0. The resulting distribution is shown in 
Fig^ We have equated the renormalization and factor- 
ization scales, fiR = fiF = fJ-, and have varied them in the 
range to/j/2 < fj, < 2to;j. There are large corrections to 
the rapidity distribution; however, they arise primarily 
from the inclusive K- factor K = ctnnlo/clo- The ra- 
pidity dependence of the NNLO X-factor is small, and is 
insignificant if normalized to the NLO distribution. This 
is not unexpected; the production of the Higgs boson at 
the LHC is strongly dominated by the soft z — > 1 limit, 
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FIG. 2: Bin-integrated rapidity distribution for LHC kinemat- 
ics, with a jet veto of < 40 GeV. We have set = mh, 
and have included the LO, NLO, and NNLO results 



which implies that the kinematics of the tree level process 
is not altered significantly by hard QCD emissions. 

In FigEI we present the rapidity distribution with a 
veto on jet activity, as described in 0|. We use the cone 
algorithm with R — 0.4, and use p v ^ to = 40 GeV; we 

therefore demand < 40 GeV. For this result we set 
nih = 150 GeV, as this cut is of particular importance 
for the H — > W + W~ decay channel. We show the LO, 
NLO, and NNLO results; we note that the LO vetoed 
rapidity distribution is identical to the inclusive LO ra- 
pidity distribution, as there is no additional radiation at 
leading order in perturbation theory. The if -factor is 
much smaller for the vetoed cross section than for the 
inclusive cross section; this can be easily seen by com- 



paring the results in FigQJand Fig|5J One reason for this 
is that the average p± of the Higgs boson increases from 
NLO to NNLO; we find < p^ LO >= 37.5 GeV, while 
< pNNLO ^_ 44 g Q e y. Since the additional partons re- 
coil against the Higgs, this is equivalent to an increase in 
the average jet p 3 ± . The increase of p 3 ± means that less 
of the cross section passes the veto, leading to a larger 
reduction in the NNLO cross-section relative to the NLO 
one. 

We have described a calculation of the fully differen- 
tial cross section for Higgs hadro-production at NNLO. 
We have utilized the method developed in pj to extract 
and cancel double real radiation singularities, and have 
extended this approach to handle processes with initial- 
state collinear singularities. Our result allows arbitrar- 
ily differential observables to be studied; as examples, we 
have presented results for the Higgs rapidity distributions 
at the LHC, both with and without a veto on jet activity. 
This approach can be easily extended to include decays 
of the Higgs, or can be applied to other 2 — > 1 processes 
of phenomcnological interest. We will provide a detailed 
description of our calculation and present a public version 
of our numerical code in a forthcoming publication. 
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